Using birth-death processes to infer tumor subpopulation structure from live-cell imaging drug screening data.

Tumor heterogeneity is a complex and widely recognized trait that poses significant challenges in developing effective cancer therapies. In particular, many tumors harbor a variety of subpopulations with distinct therapeutic response characteristics. Characterizing this heterogeneity by determining the subpopulation structure within a tumor enables more precise and successful treatment strategies. In our prior work, we developed PhenoPop, a computational framework for unravelling the drug-response subpopulation structure within a tumor from bulk high-throughput drug screening data. However, the deterministic nature of the underlying models driving PhenoPop restricts the model fit and the information it can extract from the data. As an advancement, we propose a stochastic model based on the linear birth-death process to address this limitation. Our model can formulate a dynamic variance along the horizon of the experiment so that the model uses more information from the data to provide a more robust estimation. In addition, the newly proposed model can be readily adapted to situations where the experimental data exhibits a positive time correlation. We test our model on simulated data (in silico) and experimental data (in vitro), which supports our argument about its advantages.

A comparison of mutation and amplification-driven resistance mechanisms and their impacts on tumor recurrence.

Tumor recurrence, driven by the evolution of drug resistance is a major barrier to therapeutic success in cancer. Tumor drug resistance is often caused by genetic alterations such as point mutation, which refers to the modification of a single genomic base pair, or gene amplification, which refers to the duplication of a region of DNA that contains a gene. These mechanisms typically confer varying degrees of resistance, and they tend to occur at vastly different frequencies. Here we investigate the dependence of tumor recurrence dynamics on these mechanisms of resistance, using stochastic multi-type branching process models. We derive tumor extinction probabilities and deterministic estimates for the tumor recurrence time, defined as the time when an initially drug sensitive tumor surpasses its original size after developing resistance. For models of amplification-driven and mutation-driven resistance, we prove law of large numbers results regarding the convergence of the stochastic recurrence times to their mean. Additionally, we prove sufficient and necessary conditions for a tumor to escape extinction under the gene amplification model, discuss behavior under biologically relevant parameters, and compare the recurrence time and tumor composition in the mutation and amplification models both analytically and using simulations. In comparing these mechanisms, we find that the ratio between recurrence times driven by amplification versus mutation depends linearly on the number of amplification events required to acquire the same degree of resistance as a mutation event, and we find that the relative frequency of amplification and mutation events plays a key role in determining the mechanism under which recurrence is more rapid for any specific system. In the amplification-driven resistance model, we also observe that increasing drug concentration leads to a stronger initial reduction in tumor burden, but that the eventual recurrent tumor population is less heterogeneous, more aggressive and harbors higher levels of drug-resistance.

A Comparison of Mutation and Amplification-Driven Resistance Mechanisms and Their Impacts on Tumor Recurrence.

Tumor recurrence, driven by the evolution of drug resistance is a major barrier to therapeutic success in cancer. Resistance is often caused by genetic alterations such as point mutation, which refers to the modification of a single genomic base pair, or gene amplification, which refers to the duplication of a region of DNA that contains a gene. Here we investigate the dependence of tumor recurrence dynamics on these mechanisms of resistance, using stochastic multi-type branching process models. We derive tumor extinction probabilities and deterministic estimates for the tumor recurrence time, defined as the time when an initially drug sensitive tumor surpasses its original size after developing resistance. For models of amplification-driven and mutation-driven resistance, we prove law of large numbers results regarding the convergence of the stochastic recurrence times to their mean. Additionally, we prove sufficient and necessary conditions for a tumor to escape extinction under the gene amplification model, discuss behavior under biologically relevant parameters, and compare the recurrence time and tumor composition in the mutation and amplification models both analytically and using simulations. In comparing these mechanisms, we find that the ratio between recurrence times driven by amplification vs. mutation depends linearly on the number of amplification events required to acquire the same degree of resistance as a mutation event, and we find that the relative frequency of amplification and mutation events plays a key role in determining the mechanism under which recurrence is more rapid. In the amplification-driven resistance model, we also observe that increasing drug concentration leads to a stronger initial reduction in tumor burden, but that the eventual recurrent tumor population is less heterogeneous, more aggressive, and harbors higher levels of drug-resistance.

Statistical inference of the rates of cell proliferation and phenotypic switching in cancer.

Recent evidence suggests that nongenetic (epigenetic) mechanisms play an important role at all stages of cancer evolution. In many cancers, these mechanisms have been observed to induce dynamic switching between two or more cell states, which commonly show differential responses to drug treatments. To understand how these cancers evolve over time, and how they respond to treatment, we need to understand the state-dependent rates of cell proliferation and phenotypic switching. In this work, we propose a rigorous statistical framework for estimating these parameters, using data from commonly performed cell line experiments, where phenotypes are sorted and expanded in culture. The framework explicitly models the stochastic dynamics of cell division, cell death and phenotypic switching, and it provides likelihood-based confidence intervals for the model parameters. The input data can be either the fraction of cells or the number of cells in each state at one or more time points. Through a combination of theoretical analysis and numerical simulations, we show that when cell fraction data is used, the rates of switching may be the only parameters that can be estimated accurately. On the other hand, using cell number data enables accurate estimation of the net division rate for each phenotype, and it can even enable estimation of the state-dependent rates of cell division and cell death. We conclude by applying our framework to a publicly available dataset.

Statistical inference of the rates of cell proliferation and phenotypic switching in cancer.

Recent evidence suggests that nongenetic (epigenetic) mechanisms play an important role at all stages of cancer evolution. In many cancers, these mechanisms have been observed to induce dynamic switching between two or more cell states, which commonly show differential responses to drug treatments. To understand how these cancers evolve over time, and how they respond to treatment, we need to understand the state-dependent rates of cell proliferation and phenotypic switching. In this work, we propose a rigorous statistical framework for estimating these parameters, using data from commonly performed cell line experiments, where phenotypes are sorted and expanded in culture. The framework explicitly models the stochastic dynamics of cell division, cell death and phenotypic switching, and it provides likelihood-based confidence intervals for the model parameters. The input data can be either the fraction of cells or the number of cells in each state at one or more time points. Through a combination of theoretical analysis and numerical simulations, we show that when cell fraction data is used, the rates of switching may be the only parameters that can be estimated accurately. On the other hand, using cell number data enables accurate estimation of the net division rate for each phenotype, and it can even enable estimation of the state-dependent rates of cell division and cell death. We conclude by applying our framework to a publicly available dataset.

Phenotypic deconvolution in heterogeneous cancer cell populations using drug-screening data.

Tumor heterogeneity is an important driver of treatment failure in cancer since therapies often select for drug-tolerant or drug-resistant cellular subpopulations that drive tumor growth and recurrence. Profiling the drug-response heterogeneity of tumor samples using traditional genomic deconvolution methods has yielded limited results, due in part to the imperfect mapping between genomic variation and functional characteristics. Here, we leverage mechanistic population modeling to develop a statistical framework for profiling phenotypic heterogeneity from standard drug-screen data on bulk tumor samples. This method, called PhenoPop, reliably identifies tumor subpopulations exhibiting differential drug responses and estimates their drug sensitivities and frequencies within the bulk population. We apply PhenoPop to synthetically generated cell populations, mixed cell-line experiments, and multiple myeloma patient samples and demonstrate how it can provide individualized predictions of tumor growth under candidate therapies. This methodology can also be applied to deconvolution problems in a variety of biological settings beyond cancer drug response.

Roadmap on plasticity and epigenetics in cancer.

The role of plasticity and epigenetics in shaping cancer evolution and response to therapy has taken center stage with recent technological advances including single cell sequencing. This roadmap article is focused on state-of-the-art mathematical and experimental approaches to interrogate plasticity in cancer, and addresses the following themes and questions: is there a formal overarching framework that encompasses both non-genetic plasticity and mutation-driven somatic evolution? How do we measure and model the role of the microenvironment in influencing/controlling non-genetic plasticity? How can we experimentally study non-genetic plasticity? Which mathematical techniques are required or best suited? What are the clinical and practical applications and implications of these concepts?

Exact site frequency spectra of neutrally evolving tumors: A transition between power laws reveals a signature of cell viability.

The site frequency spectrum (SFS) is a popular summary statistic of genomic data. While the SFS of a constant-sized population undergoing neutral mutations has been extensively studied in population genetics, the rapidly growing amount of cancer genomic data has attracted interest in the spectrum of an exponentially growing population. Recent theoretical results have generally dealt with special or limiting cases, such as considering only cells with an infinite line of descent, assuming deterministic tumor growth, or taking large-time or large-population limits. In this work, we derive exact expressions for the expected SFS of a cell population that evolves according to a stochastic branching process, first for cells with an infinite line of descent and then for the total population, evaluated either at a fixed time (fixed-time spectrum) or at the stochastic time at which the population reaches a certain size (fixed-size spectrum). We find that while the rate of mutation scales the SFS of the total population linearly, the rates of cell birth and cell death change the shape of the spectrum at the small-frequency end, inducing a transition between a 1/j2 power-law spectrum and a 1/j spectrum as cell viability decreases. We show that this insight can in principle be used to estimate the ratio between the rate of cell death and cell birth, as well as the mutation rate, using the site frequency spectrum alone. Although the discussion is framed in terms of tumor dynamics, our results apply to any exponentially growing population of individuals undergoing neutral mutations.

Comparison of Cell and Organoid-Level Analysis of Patient-Derived 3D Organoids to Evaluate Tumor Cell Growth Dynamics and Drug Response.

3D cell culture models have been developed to better mimic the physiological environments that exist in human diseases. As such, these models are advantageous over traditional 2D cultures for screening drug compounds. However, the practicalities of transitioning from 2D to 3D drug treatment studies pose challenges with respect to analysis methods. Patient-derived tumor organoids (PDTOs) possess unique features given their heterogeneity in size, shape, and growth patterns. A detailed assessment of the length scale at which PDTOs should be evaluated (i.e., individual cell or organoid-level analysis) has not been done to our knowledge. Therefore, using dynamic confocal live cell imaging and data analysis methods we examined tumor cell growth rates and drug response behaviors in colorectal cancer (CRC) PDTOs. High-resolution imaging of H2B-GFP-labeled organoids with DRAQ7 vital dye permitted tracking of cellular changes, such as cell birth and death events, in individual organoids. From these same images, we measured morphological features of the 3D objects, including volume, sphericity, and ellipticity. Sphericity and ellipticity were used to evaluate intra- and interpatient tumor organoid heterogeneity. We found a strong correlation between organoid live cell number and volume. Linear growth rate calculations based on volume or live cell counts were used to determine differential responses to therapeutic interventions. We showed that this approach can detect different types of drug effects (cytotoxic vs cytostatic) in PDTO cultures. Overall, our imaging-based quantification workflow results in multiple parameters that can provide patient- and drug-specific information for screening applications.

Memory in Trait Macroevolution.

The history of a trait within a lineage may influence its future evolutionary trajectory, but macroevolutionary theory of this process is not well developed. For example, consider the simplified binary trait of living in cave versus surface habitat. The longer a species has been cave dwelling, the more accumulated loss of vision, pigmentation, and defense may restrict future adaptation if the species encounters the surface environment. However, the Markov model of discrete trait evolution that is widely adopted in phylogenetics does not allow the rate of cave-to-surface transition to decrease with longer duration as a cave dweller. Here we describe three models of evolution that remove this memoryless constraint, using a renewal process to generalize beyond the typical Poisson process of discrete trait macroevolution. We then show how the two-state renewal process can be used for inference, and we investigate the potential of phylogenetic comparative data to reveal different influences of trait duration, or memory in trait evolution. We hope that such approaches may open new avenues for modeling trait evolution and for broad comparative tests of hypotheses that some traits become entrenched.

Understanding the role of phenotypic switching in cancer drug resistance.

The emergence of acquired drug resistance in cancer represents a major barrier to treatment success. While research has traditionally focused on genetic sources of resistance, recent findings suggest that cancer cells can acquire transient resistant phenotypes via epigenetic modifications and other non-genetic mechanisms. Although these resistant phenotypes are eventually relinquished by individual cells, they can temporarily 'save' the tumor from extinction and enable the emergence of more permanent resistance mechanisms. These observations have generated interest in the potential of epigenetic therapies for long-term tumor control or eradication. In this work, we develop a mathematical model to study how phenotypic switching at the single-cell level affects resistance evolution in cancer. We highlight unique features of non-genetic resistance, probe the evolutionary consequences of epigenetic drugs and explore potential therapeutic strategies. We find that even short-term epigenetic modifications and stochastic fluctuations in gene expression can drive long-term drug resistance in the absence of any bona fide resistance mechanisms. We also find that an epigenetic drug that slightly perturbs the average retention of the resistant phenotype can turn guaranteed treatment failure into guaranteed success. Lastly, we find that combining an epigenetic drug with an anti-cancer agent can significantly outperform monotherapy, and that treatment outcome is heavily affected by drug sequencing.

Glioblastoma Recurrence and the Role of O6-Methylguanine-DNA Methyltransferase Promoter Methylation.

Tumor recurrence in glioblastoma multiforme (GBM) is often attributed to acquired resistance to the standard chemotherapeutic agent, temozolomide (TMZ). Promoter methylation of the DNA repair gene MGMT (O6-methylguanine-DNA methyltransferase) has been associated with sensitivity to TMZ, whereas increased expression of MGMT has been associated with TMZ resistance. Clinical studies have observed a downward shift in MGMT methylation percentage from primary to recurrent stage tumors; however, the evolutionary processes that drive this shift and more generally the emergence and growth of TMZ-resistant tumor subpopulations are still poorly understood. Here, we develop a mathematical model, parameterized using clinical and experimental data, to investigate the role of MGMT methylation in TMZ resistance during the standard treatment regimen for GBM-surgery, chemotherapy, and radiation. We first found that the observed downward shift in MGMT promoter methylation status between detection and recurrence cannot be explained solely by evolutionary selection. Next, our model suggests that TMZ has an inhibitory effect on maintenance methylation of MGMT after cell division. Finally, incorporating this inhibitory effect, we study the optimal number of TMZ doses per adjuvant cycle for patients with GBM with high and low levels of MGMT methylation at diagnosis.

Combination therapeutics of Nilotinib and radiation in acute lymphoblastic leukemia as an effective method against drug-resistance.

Philadelphia chromosome-positive (Ph+) acute lymphoblastic leukemia (ALL) is characterized by a very poor prognosis and a high likelihood of acquired chemo-resistance. Although tyrosine kinase inhibitor (TKI) therapy has improved clinical outcome, most ALL patients relapse following treatment with TKI due to the development of resistance. We developed an in vitro model of Nilotinib-resistant Ph+ leukemia cells to investigate whether low dose radiation (LDR) in combination with TKI therapy overcome chemo-resistance. Additionally, we developed a mathematical model, parameterized by cell viability experiments under Nilotinib treatment and LDR, to explain the cellular response to combination therapy. The addition of LDR significantly reduced drug resistance both in vitro and in computational model. Decreased expression level of phosphorylated AKT suggests that the combination treatment plays an important role in overcoming resistance through the AKT pathway. Model-predicted cellular responses to the combined therapy provide good agreement with experimental results. Augmentation of LDR and Nilotinib therapy seems to be beneficial to control Ph+ leukemia resistance and the quantitative model can determine optimal dosing schedule to enhance the effectiveness of the combination therapy.

Quantifying the Dynamics of Field Cancerization in Tobacco-Related Head and Neck Cancer: A Multiscale Modeling Approach.

High rates of local recurrence in tobacco-related head and neck squamous cell carcinoma (HNSCC) are commonly attributed to unresected fields of precancerous tissue. Because they are not easily detectable at the time of surgery without additional biopsies, there is a need for noninvasive methods to predict the extent and dynamics of these fields. Here, we developed a spatial stochastic model of tobacco-related HNSCC at the tissue level and calibrated the model using a Bayesian framework and population-level incidence data from the Surveillance, Epidemiology, and End Results (SEER) registry. Probabilistic model analyses were performed to predict the field geometry at time of diagnosis, and model predictions of age-specific recurrence risks were tested against outcome data from SEER. The calibrated models predicted a strong dependence of the local field size on age at diagnosis, with a doubling of the expected field diameter between ages at diagnosis of 50 and 90 years, respectively. Similarly, the probability of harboring multiple, clonally unrelated fields at the time of diagnosis was found to increase substantially with patient age. On the basis of these findings, we hypothesized a higher recurrence risk in older than in younger patients when treated by surgery alone; we successfully tested this hypothesis using age-stratified outcome data. Further clinical studies are needed to validate the model predictions in a patient-specific setting. This work highlights the importance of spatial structure in models of epithelial carcinogenesis and suggests that patient age at diagnosis may be a critical predictor of the size and multiplicity of precancerous lesions. Cancer Res; 76(24); 7078-88. ©2016 AACR.

Optimized Treatment Schedules for Chronic Myeloid Leukemia.

Over the past decade, several targeted therapies (e.g. imatinib, dasatinib, nilotinib) have been developed to treat Chronic Myeloid Leukemia (CML). Despite an initial response to therapy, drug resistance remains a problem for some CML patients. Recent studies have shown that resistance mutations that preexist treatment can be detected in a substantial number of patients, and that this may be associated with eventual treatment failure. One proposed method to extend treatment efficacy is to use a combination of multiple targeted therapies. However, the design of such combination therapies (timing, sequence, etc.) remains an open challenge. In this work we mathematically model the dynamics of CML response to combination therapy and analyze the impact of combination treatment schedules on treatment efficacy in patients with preexisting resistance. We then propose an optimization problem to find the best schedule of multiple therapies based on the evolution of CML according to our ordinary differential equation model. This resulting optimization problem is nontrivial due to the presence of ordinary different equation constraints and integer variables. Our model also incorporates drug toxicity constraints by tracking the dynamics of patient neutrophil counts in response to therapy. We determine optimal combination strategies that maximize time until treatment failure on hypothetical patients, using parameters estimated from clinical data in the literature.

Leveraging Hypoxia-Activated Prodrugs to Prevent Drug Resistance in Solid Tumors.

Experimental studies have shown that one key factor in driving the emergence of drug resistance in solid tumors is tumor hypoxia, which leads to the formation of localized environmental niches where drug-resistant cell populations can evolve and survive. Hypoxia-activated prodrugs (HAPs) are compounds designed to penetrate to hypoxic regions of a tumor and release cytotoxic or cytostatic agents; several of these HAPs are currently in clinical trial. However, preliminary results have not shown a survival benefit in several of these trials. We hypothesize that the efficacy of treatments involving these prodrugs depends heavily on identifying the correct treatment schedule, and that mathematical modeling can be used to help design potential therapeutic strategies combining HAPs with standard therapies to achieve long-term tumor control or eradication. We develop this framework in the specific context of EGFR-driven non-small cell lung cancer, which is commonly treated with the tyrosine kinase inhibitor erlotinib. We develop a stochastic mathematical model, parametrized using clinical and experimental data, to explore a spectrum of treatment regimens combining a HAP, evofosfamide, with erlotinib. We design combination toxicity constraint models and optimize treatment strategies over the space of tolerated schedules to identify specific combination schedules that lead to optimal tumor control. We find that (i) combining these therapies delays resistance longer than any monotherapy schedule with either evofosfamide or erlotinib alone, (ii) sequentially alternating single doses of each drug leads to minimal tumor burden and maximal reduction in probability of developing resistance, and (iii) strategies minimizing the length of time after an evofosfamide dose and before erlotinib confer further benefits in reduction of tumor burden. These results provide insights into how hypoxia-activated prodrugs may be used to enhance therapeutic effectiveness in the clinic.

A high-content image-based method for quantitatively studying context-dependent cell population dynamics.

Tumor progression results from a complex interplay between cellular heterogeneity, treatment response, microenvironment and heterocellular interactions. Existing approaches to characterize this interplay suffer from an inability to distinguish between multiple cell types, often lack environmental context, and are unable to perform multiplex phenotypic profiling of cell populations. Here we present a high-throughput platform for characterizing, with single-cell resolution, the dynamic phenotypic responses (i.e. morphology changes, proliferation, apoptosis) of heterogeneous cell populations both during standard growth and in response to multiple, co-occurring selective pressures. The speed of this platform enables a thorough investigation of the impacts of diverse selective pressures including genetic alterations, therapeutic interventions, heterocellular components and microenvironmental factors. The platform has been applied to both 2D and 3D culture systems and readily distinguishes between (1) cytotoxic versus cytostatic cellular responses; and (2) changes in morphological features over time and in response to perturbation. These important features can directly influence tumor evolution and clinical outcome. Our image-based approach provides a deeper insight into the cellular dynamics and heterogeneity of tumors (or other complex systems), with reduced reagents and time, offering advantages over traditional biological assays.

An Evolutionary Approach for Identifying Driver Mutations in Colorectal Cancer.

The traditional view of cancer as a genetic disease that can successfully be treated with drugs targeting mutant onco-proteins has motivated whole-genome sequencing efforts in many human cancer types. However, only a subset of mutations found within the genomic landscape of cancer is likely to provide a fitness advantage to the cell. Distinguishing such "driver" mutations from innocuous "passenger" events is critical for prioritizing the validation of candidate mutations in disease-relevant models. We design a novel statistical index, called the Hitchhiking Index, which reflects the probability that any observed candidate gene is a passenger alteration, given the frequency of alterations in a cross-sectional cancer sample set, and apply it to a mutational data set in colorectal cancer. Our methodology is based upon a population dynamics model of mutation accumulation and selection in colorectal tissue prior to cancer initiation as well as during tumorigenesis. This methodology can be used to aid in the prioritization of candidate mutations for functional validation and contributes to the process of drug discovery.

The Impact of Microenvironmental Heterogeneity on the Evolution of Drug Resistance in Cancer Cells.

Therapeutic resistance arises as a result of evolutionary processes driven by dynamic feedback between a heterogeneous cell population and environmental selective pressures. Previous studies have suggested that mutations conferring resistance to epidermal growth factor receptor (EGFR) tyrosine kinase inhibitors (TKI) in non-small-cell lung cancer (NSCLC) cells lower the fitness of resistant cells relative to drug-sensitive cells in a drug-free environment. Here, we hypothesize that the local tumor microenvironment could influence the magnitude and directionality of the selective effect, both in the presence and absence of a drug. Using a combined experimental and computational approach, we developed a mathematical model of preexisting drug resistance describing multiple cellular compartments, each representing a specific tumor environmental niche. This model was parameterized using a novel experimental dataset derived from the HCC827 erlotinib-sensitive and -resistant NSCLC cell lines. We found that, in contrast to in the drug-free environment, resistant cells may hold a fitness advantage compared to parental cells in microenvironments deficient in oxygen and nutrients. We then utilized the model to predict the impact of drug and nutrient gradients on tumor composition and recurrence times, demonstrating that these endpoints are strongly dependent on the microenvironment. Our interdisciplinary approach provides a model system to quantitatively investigate the impact of microenvironmental effects on the evolutionary dynamics of tumor cells.

Patterns of somatically acquired amplifications and deletions in apparently normal tissues of ovarian cancer patients.

Little is understood about the occurrence of somatic genomic alterations in normal tissues and their significance in the context of disease. Here, we identified potential somatic copy number alterations (pSCNAs) in apparently normal ovarian tissue and peripheral blood of 423 ovarian cancer patients. There were, on average, two to four pSCNAs per sample detectable at a tissue-level resolution, although some individuals had orders of magnitude more. Accordingly, we estimated the lower bound of the rate of pSCNAs per cell division. Older individuals and BRCA mutation carriers had more pSCNAs than others. pSCNAs significantly overlapped with Alu and G-quadruplexes, and the affected genes were enriched for signaling and regulation. Some of the amplification/deletion hotspots in pan-cancer genomes were hot spots of pSCNAs in normal tissues as well, suggesting that those regions might be inherently unstable. Prevalence of pSCNA in peripheral blood predicted survival, implying that mutations in normal tissues might have consequences for cancer patients.